Mathematical Issues in a Fully-Constrained Formulation of Einstein Equations 



Isabel Cordero-Carrion/'El Jose Maria Ibanez/'QEric Gourgoulhon, 2, [|| Jose Luis Jaramillo, 3 ' 2 'Hand Jerome Novak 2 '0 

1 Departamento de Astronomia y Astrofisica, Universidad de Valencia, Valencia, Spain 
2 Laboratoire Univers et Theories (LUTH), Observatoire de Paris, CNRS, 
Universite Paris Diderot, Place Jules Janssen, 92190 Meudon, France 
"' Instituto de Astrofisica de Andalucia, CSIC, Apartado Postal 3004, Granada 18080, Spain 

(Dated: 18 February 2008) 

Bonazzola, Gourgoulhon, Grandclement, and Novak [Phys. Rev. D 70, 104007 (2004)] proposed 
a new formulation for 3+1 numerical relativity. Einstein equations result, according to that for- 
malism, in a coupled elliptic-hyperbolic system. We have carried out a preliminary analysis of the 
mathematical structure of that system, in particular focusing on the equations governing the evolu- 
tion for the deviation of a conformal metric from a flat fiducial one. The choice of a Dirac's gauge for 
the spatial coordinates guarantees the mathematical characterization of that system as a (strongly) 
hyperbolic system of conservation laws. In the presence of boundaries, this characterization also 
^ ' depends on the boundary conditions for the shift vector in the elliptic subsystem. This interplay 

\ between the hyperbolic and elliptic parts of the complete evolution system is used to assess the 

prescription of inner boundary conditions for the hyperbolic part when using an excision approach 
to black hole spacetime evolutions. 
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D" 1 I. A FULLY-CONSTRAINED EVOLUTION SCHEME 

U . 

A second-order fully-constrained evolution formalism for the Einstein equations has been proposed in Ref. [l£ 
This evolution scheme, that will be referred in the following as Fully-Constrained Formulation (FCF), is based on 
a conformal 3+1 formulation of General Relativity and makes use of an elliptic condition for the choice of spatial 
coordinates, a generalized Dirac gauge, and a maximal condition for the slicing. The enforcement of the constraints 
along the evolution together with the elliptic nature of the employed gauge conditions, translates the FCF formalism 
into a mixed elliptic-hyperbolic Partial Differential Equations (PDE) system, consisting in five quasi-linear elliptic 
equations coupled with a tensorial second-order in time and in space evolution equation for the conformal metric. In 
this article, we aim at gaining insight on some mathematical issues associated with this PDE system and, in particular, 
assessing the hyperbolicity of the tensorial evolution part. A good understanding of the mathematical structure of 
the system will be crucial in the context of full 3D numerical relativity simulations, since the choice of state-of-the-art 
numerical tools will be adapted to the specific structures of the whole system governing the evolution of matter fields in 
• • ■ a dynamical space-time: spectral methods for the elliptic subsystem [32|, and modern high-resolution shock-capturing 
. £h ' techniques for the hyperbolic part [H, H3| ■ The implementation of the scheme in [l8[ will naturally extend previous 
works — following the Conformal Flatness Condition (CFC) approach of Isenberg- Wilson-Mathews [H, H3] — devoted 
5_i ' to the study of some relevant astrophysical sources of gravitational radiation [2l], [22|, [23|, [24[ • 

A. Gauge reduction, PDE evolution systems and well-posedness 

The gauge character of General Relativity (GR) strongly conditions any attempt of finding a solution by solving a 
Partial Differential Equations (PDE) problem. In its standard formulation through the Einstein equation 

Rtiu- -Rg llv = ^T l j, v , (1) 
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solutions are given in terms of spacetime geometries (A4,g^), i.e. classes of Lorentzian metrics g^ v equivalent under 
diffcomorphisms of M. , rather than by specific 4-metrics in some particular coordinate system. As a consequence of 
this, any attempt to cast ([1]) as a standard PDE system necessarily must go through a gauge reduction process. This 
fixing of the gauge involves four different (differential) systems: i) the reduced system, whose solution provides the 
metric in a given coordinate system, ii) the constraint system, consequence of the gauge character of the theory and 
that characterizes the solution manifold, iii) the gauge system, which fixes the coordinate chart and permits to write 
the reduced system as a standard PDE problem, and iv) the subsidiary system, guaranteeing the overall consistency 
along the evolution and, in particular, between the reduced and gauge systems. The mathematical consistency of 
the evolution formalism involves two aspects. First, one must assess the analytic well-posedness of the PDE system 
that is actually solved during the evolution, that we will refer to in the following as the evolution PDE system, that 
includes the reduced system but possibly other additional PDEs. Second, one must guarantee the fulfillment of the 
subsidiary system during the evolution. 

As in other evolution formalisms based on the Initial Value problem for the Einstein equation [2j, the constrained 
system in the FCF scheme follows from the Gauss-Codazzi-Ricci conditions 

- K i: jK tj + K 2 = 16np 
Dj (K ij - j ij K) = 8n J 1 , (2) 

i.e. the Hamiltonian and momentum constraints in the 3+1 formulation (p is the energy density and J 1 the current 
vector) which are elliptic in nature. The currently most successful numerical evolution formalisms are free schemes in 
which the constraint system is not enforced during the evolution. This is the case of certain generalized harmonic 
formalisms [!, 0| and the 3+1 BSSN (from Baumgarte, Shapiro, Shibata and Nakamura; see references @,1) used 
in recent binary black hole breakthroughs @, H, 0, G3 and in fully 3D evolution of binary neutron stars (see e.g. 
[HI). In these free schemes, the corresponding evolution PDE system is formed by the respective reduced systems 
together with some additional evolution equations to fix the harmonic gauge sources, in the case generalized harmonic 
schemes, or the lapse function and shift vector, in the BSSN case. No elliptic equation is solved during the evolution 
and standard hyp erbolic techniques can in principle be used to assess the well-posedness of the evolution system (cf. 
in this sense [12| for the case of the BSSN system). In contrast, the FCF here discussed actually incorporates the 
constraints to the evolution PDE system. Moreover, the use of the above-mentioned elliptic gauge conditions adds 
additional elliptic equations during the evolution. The resulting FCF scheme presents some interesting properties as 
compared with free evolution schemes. Apart from the absence of constraint violations (an issue under control in 
current BSSN and generalized harmonic formulations), we can highlight the following features (cf. [l8[ for a more 
complete discussion): first, the FCF naturally generalizes (as commented above) the successful scheme employed 
in the CFC approximation to General Relativity; second, it permits to read the gravitational waveforms directly 
from the metric components; third, the scheme can be straightforwardly adapted to the extraction of gravitational 
radiation at null infinity by making use of hyperboloidal 3-slices implemented by means of a constant mean curvature 
elliptic gauge condition; and fourth, it provides a well-suited framework for the formulation of realistic (approximate) 
prescriptions in the construction of quasi-stationary astrophysically configurations [HI]. However, the well-posedness 
analysis of such a mixed elliptic-hyperbolic system can be a formidable problem, since part of the dynamics related to 
the characteristic fields in the hyperbolic part is encoded in fields obtained only once the elliptic part is solved. Even 
though analyses of such systems exist in the GR literature (see e.g. Refs.fl^. fl5L [l6| and particularly Ref. they 
deal with free evolution systems, in which the elliptic part follows only from the gauge conditions. The well-posedness 
analysis of the complete elliptic-hyperbolic system in the FCF scheme, which in addition includes the constraints, 
is beyond the scope of this work and we will mainly focus on the hyperbolicity analysis of the tensorial evolution 
equation. Before referring to the additional issues related to the subsidiary system, we must provide some details 
about the FCF formalism. 



B. Brief review of the FCF scheme 

Following Ref. [Ts| . we consider a standard 3+1 decomposition of an asymptotically flat spacetime {M^g^v) in 
terms of a foliation by spacelike hypersurfaces (St). We denote the unit timelike normal vector to the spacelike slice 
S t by n^, the spatial 3-metric by 7^, i.e. 7 M „ = g^ v + n fJi n v , and adopt the following sign convention for the extrinsic 
curvature: K^ u = —^Cnj^. The evolution vector t^ = (9t) M is decomposed in terms of the lapse function TV and 
the shift vector (3^, as t^ = Nn^ + (3^. 

Under this 3+1 decomposition, Einstein equation (TT]) splits into the 3+1 constraints in (|2|) and a set of evolution 
equations for the extrinsic curvature that, together with the kinematical relation defining the extrinsic curvature, 
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constitute the 3+1 evolution equations 

(d t - Cp) jij = -2NK t] 

(fit - C ) Kn = —D t DjN + N { ^ + KKij - 2K k K k] + 4^ [(S - £) 7ij - 2S l3 ] } . (3) 

This is a first-order in time and second-order in space evolution system for ( 7i j , K 1J ) . 

The first specific element in the FCF scheme is the introduction of a time independent fiducial flat metric fy, 
which satisfies Ctfij — dtfij = 0. This rigid structure is chosen to coincide with 7^ at spatial infinity, capturing its 
asymptotic Euclidean character, and permits to work with tensor quantities rather than with tensor densities. We 
will denote by T>i the Levi-Civita connection associated with . 

a. Conformal decomposition. As a step forward in the reduction process to the PDE system in the present FCF, 
we perform a conformal decomposition of the 3+1 fields: 

7«=* 4 7« ^ = * 4 i« + \K^ , (4) 

where K — 7 y i<Q J -, the representative 7^ of the conformal class of the 3-metric is chosen to satisfy the unimodular 
condition det( 7 y ■) = det(fij), and the traceless part A lJ of the extrinsic curvature is decomposed as 



A ij = ^ + &p - \b k p k fi + d t f^ , 



(5) 



with Di the Levi-Civita connection associated with 7 y . Finally, in the following we will denote by h lJ the deviation 
of the conformal metric from the flat fiducial metric, i.e. 

h ij := f j - f ij . (6) 

Using these conformal decompositions of 7^ and K % \ the 3+1 constraints @ and evolution system ^ can be 
expressed in terms of the basic variables h % \ "J, N, (3 l , K. Before giving more explicit expressions, let us remove the 
gauge freedom. 

b. Gauge system. Following the prescriptions in [l8j |. namely maximal slicing and the so-called generalized Dirac 
gauge, we choose 

K = 0, iP :=V k ^ = 0, (7) 

These g aug e conditions fix the coordinates, even in the initial slice, up to boundary terms (see e.g. sections 9.3. and 
9.4. in [23). These two relations define the gauge system in the FCF scheme. Since the gauge system is meant to 
hold at all times, the following conditions must also be satisfied 

k = o, d t (v k f i ) = 0. (8) 

The FCF scheme actually enforces the first of these conditions, K — 0, during the evolution. Taking the trace in the 
second equation in ([3]), and using the Hamiltonian constraint that is also enforced during the evolution (see below), 
an elliptic equation for the lapse follows 

D k D k N + 2D k In * D k N = S N [N, /J\ (9) 

c. Main or reduced system. In the FCF scheme in Ref. 18] the reduced system is a second-order in time and 
second-order in space evolution system for the deviation tensor h lJ . This is obtained by: i) combining equations in 
([3]) into a single second-order in time equation; ii) inserting in it the conformal decompositions ^ and ([5]), and iii) 
imposing the gauges ([7]). The resulting expression is formally written as (see next section for a detailed account): 

~ir - ^ kl vM^ - 2C^ + CeCpH<* = 3%, (10) 

where the source S 1 ^ does not contain second derivatives of . Use of the Dirac gauge results in the wave-like form 
of this equation, since it eliminates certain second derivatives of the type V 1 !)^^ coming from the expression of the 
Ricci tensor. 
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d. Constrained system. The Hamiltonian constraint in ([2]) can be written as an elliptic equation for the conformal 
factor 

D k D k V~ 3 -fy = S*[tf,jV,/? i ) 7 ij ]. (11) 

Again 5* ['J, N, f} % , 7^] represents a non-linear source. Momentum constraint poses a more subtle issue. In Ref. [l8| 
an elliptic equation for the shift vector is deduced using both the momentum constraint and the preservation in time 
of the Dirac gauge (second relation in ©): 

D k D k p + ^D l D k (3 k + 3 R l k (3 k = S£[¥, N, (3\ %} (12) 

An equation for the shift could be derived from the momentum constraint alone, but the coupling to the tensorial 
equation (fl"T)|) would become more complicated due to the presence of a mixed time-space second-order derivative of 
h % i . This term is eliminated by the use of a Dirac, or a similar, gauge. 

Alternatively, an elliptic equation for the shift can be drawn from the preservation of the Dirac gau ge a lone, 



renouncing, therefore, to the fully-constrained character of the scheme — e.g. this is the strategy in Ref. jl7|, but 
using a spatial harmonic gauge condition instead of the Dirac one. At the end of the day, the choice (fT2")) in the FCF 
scheme provides an elliptic equation for the shift that enforces the momentum constraint, as long as the Dirac gauge 
is satisfied. 

e. FCF evolution PDF system. The mixed elliptic-hyperbolic PDE system that evolves some initial data given 
on an Cauchy slice is formed by: a) Eqs. ©, (fTI]) and (TT2")) . the elliptic part, and b) Eq. (|10p . the wave-like tensorial 
equation. As we have pointed out, we will not consider here the well-posedness analysis of the whole system. To 
give an idea of the involved difficulties, we note that the elliptic part is very similar to the Extended Conformal Thin 
Sandwich (XCTS) employed in the construction of initial data, though here it is solved all along the evolution. 

Even the restriction to the elliptic subsystem represents a very hard problem, as it is illustrated by the lack of the 
existence results for the XCTS system and the preliminary numerical [28| (see also [29j]) and analytical [3(| [M| results 
pointing towards a generic non-uniqueness of the elliptic system. For these reasons, we will focus on the study of 
the hyperbolicity of the tensorial evolution equation (|10[) , understanding this as a necessary condition for the overall 
well-posedness. 

/. Subsidiary system. The resolution of the PDE evolution system only guarantees the consistency between the 
reduced and gauge systems as far as the slicing condition is regarded, since equation @ for the lapse is indeed 
enforced. This is in principle not the case for the Dirac gauge. More dramatically, if the Dirac gauge is actually 
not satisfied, the FCF scheme is not really fully-constrained, since in that situation Eq. (|T2|) no longer enforces the 
momentum constraint. A control of the evolution of the Dirac gauge is therefore crucial in the scheme. A wave-like 
equation for T) k h k% can be obtained by taking the divergence of the tensorial Eq. (|10[) . The vanishing of T> k h kl in the 
evolution would then follow from the initial conditions V k h k% — and d t {V k h kl = 0) = imposed in the construction 

of the initial data, and the satisfaction of Eq. (91) in Ref. 18] for f3 l . The latter can be considered as the subsidiary 
system in the FCF scheme. 



C. Specific objectives and organization 

Though the wave character of Eq. (fTU|) essentially guarantees its hyperbolicity, we aim here at developing a more 
detailed analysis. This is motivated by the need of controlling the characteristics in initial boundary problems and also 
when trying to make use of first-order techniques employed in matter evolutions. Our main specific goal in this article 
is the development of a hyperbolicity analysis of a first-order version of the evolution part in the FCF formalism, 
where N, VP and (3 % are considered as fixed parameters. In particular, we aim at obtaining explicit expressions for 
the characteristic fields and speeds. As pointed out above, this point represents a fundamental ingredient in the 
study of the appropriate boundary conditions if boundaries are present in the integration domain. This constitutes 
only a preliminary study of the well-posedness of the evolution system since no stability analysis whatsoever will be 
considered. Certainly further analysis is required. However, in the absence of a full treatment and being ultimately 
motivated by practical numerical implementations needs, the level of rigor and completeness in this article is adapted 
to the achievement of limited but concrete results. 

On behalf of self-consistency, and in spite of the lack of a fully rigorous treatment of the FCF subsidiary system, 
we also aim at discussing certain (numerical) algorithms devised to guarantee the fulfillment of the Dirac gauge along 
the evolution. Though this is not the substitute of a formal proof it provides, on the one hand, support for the 
coherence among the reduced, gauge and constrained systems. On the other hand, and more importantly from a 



5 



practical point of view, the implementation of the FCF scheme is then guaranteed to be fully-constrained, even in 
numerical implementations where errors can occur even if analytic well-posedness has been established. 

The article is organized as follows. Section [TT] presents first-order formulation of the FCF scheme, more concretely 
of its reduced system. In section IIIII the characteristic structure of the reduced system is analyzed, with a brief 
application to inner boundaries in excised black hole spacetime evolutions. Section IIVI discusses the possibility of 
writing the first-order reduced FCF system as a system of conservation laws, by making explicit use of the Dirac 
gauge. In section [V] two different manners of enforcing the Dirac gauge in the evolution are introduced, providing 
key support for overall consistency and guaranteeing the fully-constrained character of the scheme. Finally section 
IVII concludes with a discussion of the results. 



II. FIRST-ORDER REDUCTION OF THE REDUCED SYSTEM IN THE FCF 



Equations governing the evolution of h v in the FCF are 

d 2 h lj N 2 



dt 2 i/r 



—^V k V^ - 2C a 



dt 

Cp 1 h ij 



N 



6 V k Q (V l h jk + V ] h lk - V k h ij ) 



dt 



Cf3 I InN 



d_ 

dt 



-V k (3 k W + (LP)' 1 



0_ 

dt 



C 3 ) V k f3 k - H {p k S3 k ) 2 



h'J 



-2Nip- 4 Z lj 



-{2N) Z 
-2N^~ 6 



j kl A ik A^ - 4tt hp 4 S ij - -Sf j 
j ik j jl V k ViQ + X - (h ik V t h l i 



+h jk V k h il - h kl V k W) ViQ - -f^ kl V k V^ 



(13) 



where S ,y and S are, respectively, the spatial components of the stress tensor S a /3 :— JaJpT^, associated with the 
matter energy-momentum tensor T„ v , and its trace. {Lfl) 13 is the conformal Killing operator associated with the flat 
metric acting on the vector field j3 l : 



(Lpf := v v + n- r - -n, f\r 



(14) 



and the auxiliary quantities Q and Z 11 are 



Q := N^ 2 , 



(15) 



= N RV + 8ip~ 2 (f k V k ip) (^V^) 
+A4,- 1 (f k V k iP) (f%N) 
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--^- 1 V k ^(t l 'DhN)f j . (16) 



The symmetric tensor Rj is defined by 



K j := \ [--Dih ik V k W l - %a mn V m h lk V n V l 

+\i lk i 3l v k h mn va mn , (17) 



and the scalar i?* is 



R. := -^ kl V k h mn Va mn - ^j kl V k h mn V n j mn . (18) 

Let us write Eqs. (|13p as a first-order system, by introducing the following auxiliary variables: 

dh ij 



dt 



(19) 



J k 

With these new variables the system for h lJ can be cast into 

duV N 2 



w l ':=V k W. (20) 



— - ^ kl V k wf - 2(3 k V k u 13 + f3 k (3 l V k w\ 3 



= <F (/3 fe , TV, ^, d,(3 k , d,N, drf, , (21) 

where (j? 3 are source terms which do not contain partial derivatives of u %3 or w£ . From definition (|20[) we obtain 

= , (22) 

where we have taken into account that c't/ 1 -' =0. In terms of the above new auxiliary variables, the system of Eqs. 
P9|21|22[) . can be written as: 



~dt 



+ A%v = g (>, N, i>, d^ k ,d^N, d^, h l3 ,u l3 lW fj , (23) 



where the vector v is: 



/ (h ij 
[u ij 

(«#) 



(24) 



and the source g is 



(0 k ,N,1>,dp0 k ,drN,d^,hV^,wit} = . (25) 



(0) 

In these equations, v and g are vectors of dimension 30, as it results from the symmetry properties of h lJ , u 13 , and wV . 
Let us remind that, besides the above symmetry properties, the following algebraic constraints have to be satisfied: i) 
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det 7ij = det ; and vjf = 0, which is equivalent to Dirac's gauge. In order to write the matrices of the system 
in a simple way, the following auxiliary quantities are defined: 



q 13 := -N 2 ip- 4 f j , 

Qi ._ ( q li q 2i q 3i J ^ 

-si' 
-si 



Then, the explicit form of the matrices A' are: 

/ 06x6 



A' = 



J6x24 



-2/?'l 8 



03 



x5 



J3x4 



3 X3 



J15 



J12 



J3x2 



Q l 
Q z 







18x18 



(26) 
(27) 

(28) 



(29) 



III. CHARACTERISTIC STRUCTURE OF THE REDUCED SYSTEM 



Let us present here a preliminary analysis of the mathematical structure of system (|23[) . 

First, we give the explicit expressions of the characteristic speeds in terms of the functions ip, N, {3 l and 7 U . 
Lemma 1: Let us consider the evolution vector d t , whose components are = (1,0,0,0), and a generic spacelike 
covector of components £ Q = (0,Q) orthogonal to the evolution vector. The associated eigenvalue problem (see, e.g., 
ref. f!/j: 

[A'Ci - AI] X A = 0, (30) 
where X denotes the eigenvalue and ~X.\ the corresponding eigenvector, has the following solution: 

A = 0, 

= -^±N(ec^ 1/2 , (31) 

where Xq has multiplicity 18, and each X±^ has multiplicity 6. 

Imposing Dirac's gauge in (|7|) indeed guarantees the real character of the eigenvalues corresponding to matrices A 1 , 
and therefore the hyperbolicity of the evolution system. Even though this is not a prerogative of the Dirac gauge, 
other prescriptions for H l in condition ([7]) lead to a more complicated structure of the resulting sources. As mentioned 
after Eq. (|12p. a more important point is the fact that other choices of H l will generally introduce time derivatives 
of in the elliptic subsystem, complicating further the complete PDE system. Of course, if no gauge is imposed at 
all, one can check that the A' matrices admit complex eigenvalues. This reflects the property that Einstein equations 
by themselves do not have a definite type, without the specification of a gauge. We conclude that when imposing 
Dirac's gauge the eigenvalues of the linear combination A l Q are real: 
Lemma 2: Dirac 's gauge is a sufficient condition for the hyperbolicity of system (|l?ffp . 

In the above eigenvalue problem, the first 6 eigenvectors, with eigenvalue and associated with the h 13 components 
of v in (f2"l)) . completely decouple from the other eigenvectors. Therefore, the rest of eigenvectors can be studied 
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independently. For the sake of clarity in the notation, let us define some auxiliary quantities before writing the matrix 
of (right-)eigenvectors: 



Ci 



C 2 := 



-Qq 12 ~( t q l3 \ 
Qq 11 
&q a J 



VCa 



(32) 



The matrix of (right) eigenvectors, R,^, associated with the eigenvalue problem described in the above Lemma 1 



is: 





( h 


06x24 






06x12 






R<« = 


024x6 


Ci 







c 2 

c 2 





c 2 




c 2 




V 












c 2 





c 2 



(33) 



If the determinant of this matrix vanishes, the set of eigenvalues is not complete. This happens in the following cases: 



Case 1: \ { P = \ {0 . Since 



(34) 



and does not vanish (£* is a spatial vector different from zero) non-completeness only occurs if the lapse N 
vanishes. 

Case 2: (iQq 13 = 0. From the definition of q v , it follows 



(35) 



One can see that the previous equality depends only on the direction of the vector (i.e. = 1). From now 
up to the end of the study of the different cases, the vector Q will be considered to be unitary. So ([35]) loads to: 



QQ (FP - WV- 4 7«) = & (bp) 2 = N 2 



(36) 



Decomposing ft 1 into components parallel and normal to C\ we write (3 l = (/?") C, 1 + , where (/?") = Q(3 l 

and Ci (Z^)' = 0. From flU), we conclude: 



GC,V J ' = 0^ (/?")' 



TV 2 



(37) 



Note that this case is independent of the choice of £*, since it corresponds to (/?")* (^") ., i.e. |C*A| 2 - Therefore, 
non-completeness occurs if |/?"| = N. 

Case 3: Qq lJ = 0, Vj = 1, 2, 3. This is a stronger case than the previous one. Again from the definition of g u , 
we have: 



& - N 2 ^-^) = & (C^) ft = N 2 (i 
From this, and the decomposition (3 l — (/3") ( l + it follows: 

Qq 13 =0<& (13^ =0 and ^ 2 = N 2 
This is just a stronger version of the second case above. 



(38) 



(39) 
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As a consequence of the above analysis we can set up the following lemma. 

Lemma 3: The (right-) eigenvectors associated with the matrix A l Q define a complete system iff i) the lapse N does 
not vanish, and ii) the projection of the evolution vector onto the plane spanned by andQ^, i.e. = Nn^ + & / 3£ M , 

is non-null, i.e. (/?") ^ N 2 . 

In the eigenvalue problem (l30|) . Q stands for an arbitrary spatial vector. In particular, we can always choose C — (3 l . 
In that case, the degeneracy condition in cases 2 and 3 above reduces to (3 % j3i = N 2 . This happens if the vector f 
becomes null. Moreover, if the vector t^ is spacelike then we are in case 2, since then there exists a vector (in fact, 
a cone obtained by the rotation of the non-vanishing j3 l by an appropriate angle) such that the projection of P l onto 

that C, refered to as (/?H)\ satisfies (fS^f (^). = (Qft) 2 = N 2 . We conclude: 

Proposition 1 : The system H23]) is strongly hyperbolic if 'f is timelike, i.e. if N ^ and N 2 — (3 l > 0. 



In some particular cases, degeneracy in the eigenvalues can occur. In particular, it could happen that one of the 
eigenvalues A + or A_ coincides with A . These degeneracies can appear where: 

aL° = O (/3%,) 2 = N 2 (C%) • (40) 

Again, one can consider to be unitary. Hence, either A+ or A_ vanishes when (/?") 2 = N 2 . As seen in (f36|) , in this 
case the system of eigenvectors is incomplete. 

Another relevant property is the following: 
Proposition 2: All the characteristic fields associated with the eigenvalue problem \3Cfy are linearly degenerate, i.e., 
they satisfy the following condition: 

DA p (v) ■ r p (v) = 0, (41) 

where r p is the eigenvector associated to the eigenvalue X p , and the operator D is defined in the space of the variables 
of the system. 

This shows the good behaviour of the Dirac gauge since, in the language from fluid dynamics, it means that no shocks 
can be propagated along these curves, in particular gauge shocks. Hence, if there were discontinuities, they have to 
be contact discontinuities. 

Regarding the characteristics speeds A± we have: 
Corollary 1: The non-zero eigenvalues associated with Q l correspond to the coordinate velocity of light. 
This feature, which is an expected result, can be shown by considering a unitary Q l and a curve whose spatial part 

,/.,•' 

' £*. Using the 3+1 expression of the metric, the vanishing of the line element 



dx % 

points in the C direction: — — 

dt 



dt 



of the curve, where the component of j3 l in the direction is considered, is imposed. It follows, using the expression 

dx 



for A^ in (ED that A (c) 



dt 



A. Application to inner boundary conditions 

The explicit expressions (|3"Tj) for the characteristic speeds are specially useful in the assessment of the boundary 
conditions to be imposed on a given border. We illustrate this by considering inner boundaries in the context of excised 
black hole spacetimes. Before doing so, let us underline that the FCF can be employed in combination with any of the 
standard techniques dealing with the black hole singularity in numerical evolutions of black hole spacetimes, namely 
excision, punctures or stuffed black holes. However, the excision technique is favoured if (the elliptic subsystem of) the 
FCF is implemented by means of spectral methods. Focusing on the excision approach, let us denote by <S f the inner 
sphere employed as inner boundary at a given spacelike slice £(, and by TL the worldtube hypersurface generated 
along the evolution by piling up the different St- A natural expectation is that no inner boundary conditions should 
be prescribed for radiation fields on inner supcrluminal (growing) inner boundaries. This would avoid the need to 
incorporate boundary conditions in the well-posedness analysis of the associated initial boundary value problem. From 
this reason, spacelike inner hypersurfaces TL are good candidates for inner boundary conditions. However, this general 
idea must be assessed in the context of every specific evolution scheme. In our particular case, we must check that 
characteristic speeds (l3~Tj) arc outgoing (with respect to the integration domain). The tangent vector h^ to TL which 
is normal to each St , and transports St into St+st , can be written as 



hP = Nn^ + hss^, 



(42) 
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where s M is the normal vector to St, lying on T, t and pointing toward spatial infinity. Then, since the norm of is 
given by h^h^ — —N 2 + h 2 s , it follows that H is spacelike as long as b > N. Choosing a coordinate system adapted 
to H, i.e. where all the spheres St stay at the same coordinate position — say r = const = r — it follows that 
h s = (3 l Si = (3- 1 . In this case, H is spacelike as long as (3 1 - > N. Evaluating expression ([3"T]) for £' = s l , it follows 

\ { ± ] = -p 1 - ± N (43) 

From this it follows that: 

Corollary 2: For a coordinate system adapted to a spacelike inner worldtube TC, where /3 > N, no ingoing radiative 
modes flow into the integration domain S t at the excision surface. 

Under these conditions no inner boundary conditions whatsoever must be prescribed for the hyperbolic part. Of 
course, it is not obvious how to choose dynamically an inner boundary Ti. that is guaranteed to be spacelike during 
the evolution. A proposal in this line has been presented in [371 ] in the context of the dynamical trapping horizon 
framework (see e.g. Ref. [36|). Quasi- local approaches to black hole horizons aim at modeling the boundary of a 
black hole region as world-tubes of apparent horizons (St). Dynamical horizons provide a geometric prescription for 
7i that is guaranteed to be spacelike, as long as the black hole is dynamical, and remain inside the event horion, 
if cosmic censorship holds. The corresponding geometric dynamical horizon characterization is enforced as an inner 
boundary condition on the the elliptic part of the FCF, in particular on the shift equation (fl"2")) . This shows the 
key interplay between elliptic and hyperbolic modes in the coupled fully-constrained PDE evolution system. Note 
however that, according to Proposition 1, the hyperbolic evolution system ceases to be strongly hyperbolic. In fact, 
the evolution vector < M , tangent to H. in the adapted coordinate system, becomes spacelike in a finite region. This can 
be bypassed by adopting a coordinate system in which the coordinate radii of the S t slices grow in time: r = r(t) ^ 0, 
where r(t) is appropriately chosen. In this case, h s — j3 holds no longer, and this relation is rather substituted by 
/3 1 - = h s — [r(t) — r a }. This condition is again under control through the appropriate boundary condition on the elliptic 
equation for (3 l . Note that in this case the characteristics are still outgoing from the integration domain though, in 
this case with a coordinate growing excision sphere, this feature is no longer characterized by the negativity of the 

( s) ( s) 

characteristics speeds \± . The outgoing character is guaranteed by the characterization of \± as the coordinate 
velocity of light in Corollary 1, together with the spacelike character of Tt. 

IV. DIRAC GAUGE AND SYSTEM OF CONSERVATION LAWS 

A hyperbolic system of conservation laws, without sources, is: 

<9 t u + D, ; F(u) =0. (44) 

In this system we can identify the set of unknowns, i.e., the vector of conserved quantities u, and their corresponding 
fluxes f(u). 

The choice of Dirac's gauge allows us to find the following set of I vector fluxes f' (I = 1, 2, 3), of dimension 30: 



/ (0 6 ) \ 



(45) 



in terms of which system (|23p can be rewritten as a hyperbolic system of conservation laws (with sources). The 
Jacobian matrices associated to the fluxes f', (A*)' are: 
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J6x30 





-W l h 


Q l 
Q l 


0l8x6 


-s 


03x5 


Oi8xi8 


Olfi 


-5 l 


03x4 


Ol2 


-s l 


3 X3 


Og 


-5 l 


03x2 


o 6 


-s l 


3 


03 





where 



(46) 



E' 



w[ J S[ w\[8 2) w^S^ w l 2 J 8 2 w^5 l 3) w l 3 J S 3 

E 12,l 
E 13,l 
J?22,l 
£23,1 

V E 33 ' 1 J 



(47) 



and the parentheses in the subindices represent a symmetric sum, e.g., Wp^j = W\ 5 2 + w 2 S[. 

These matrices have the same eigenvalues as the matrices A'. The corresponding eigenvectors are different but 
they keep the same fundamental properties as the ones associated to the matrices A 1 , namely they define a complete 
system. Hence, the following lemma is in order: 

Proposition 3: Taking advantage of Dirac 's gauge, it is possible to convert the hyperbolic part of the coupled elliptic- 
hyperbolic system of the FCF formalism, into a (strongly) hyperbolic system of conservation laws (with sources). 



V. PRESERVATION OF THE DIRAC GAUGE IN THE EVOLUTION: THE DIRAC SYSTEM 

The importance of the enforcement of the Dirac gauge during the evolution in time has already been stressed in 
the introduction. In this section we give a brief description of some numerical algorithms that can be used to fulfill 
the Dirac gauge, when solving the reduced system (jlOp . In particular, we do not intend to provide a formal proof 
of the consistency of the method. Because of the unimodularity of the conformal metric 7^ , the symmetric tensor 
h lJ has only five degrees of freedom. For simplicity, here we shall illustrate the scheme by considering the case where 
the trace h = fijh lJ = 0. The unimodular condition would be satisfied by an iteration on the value of the trace, as 
described in We consider the particular case of spherical polar coordinate system (r,9,(p), and note by A the 
flat Laplace operator, i.e. 

d 2 2d 1 

A :=V{D l = — + -— + — A ev , (48) 
or r or r z 

where Ag v involves only angular derivatives. Thus, the problem to be solved can be written as a wave equation with 
constraints 

£- A )^ = Sl3 > ( 49 ) 

V 3 h lj = , (50) 
h = ; (51) 

where the source iS u gathers all the other terms of Eqs. (fH?)) , including the shift terms in the differential operator. The 
structure of the differential operator in the left-hand side is here simplified with respect to the full evolution one of 
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Sec.ini in order to focus on the propagation aspects, which are already contained in the simple wave operator. The full 
evolution operator can also be handled with a similar technique, but involving more technical justifications. The sys- 
tem (j49j) - (f5Tj) can be seen as the evolution of two scalar fields, two dynamical degrees of freedom, from which one recov- 
ers the full tensor h lJ using the trace and divergence-free conditions. To gain insight, it is helpful to decompose the ten- 
sor on a basis of Mathews-Zerilli [38l.l39j tensorial spherical harmonics. We use the basis of six families of pure-spin ten- 
sor harmonics as referred to by Thorne j^, with the same notations: T L °> em , T El ' em , T Sl ' £m , T E2 > em , T B2 < em , T T °' em . 
If we note the coefficients of h lJ in this basis (c L °^ m , c El > lm , c Bl ' em , c E2 ' im , c B2 ' em , c T "- em ^ , we can define for any rank 2 
symmetric tensor the following six scalar fields: 

L := Y, cL0 ' emY *™ = hrr > 
r, := £ c E ^Y im , 

l>l,m 

p := Y, cBl/mY t™, 
e>i,m 

W : = E cE2 ' emY im, 
t>2,m 

X := ^' lm Y^ 

t>2,m 

T := £> T °'*"r /m > ( 52 ) 

where Yi m (9, ip) are the scalar spherical harmonics, which are eigenfunctions of the angular Laplace operator Ag v Ye m = 
—£(£ + l)Yf m . Note that there is a one-to-one relation between the six components of h lJ and these six scalar fields. 
The trace condition ([51) simply turns into To + h rr = 0, therefore we shall replace To with — h rr in all forthcoming 
expressions. The divergence-free conditions (|50j) turn into: 

g h rr 3h rr i 

-&e v V = 0, (53) 



Or 



r r 



b + h + {Al , e + 2) ™Jf = l) , (54) 

or r r Zr 

^ + ^ + (A e , + 2)^ = 0; (55) 
or r r 



where all the angular derivatives are expressed in terms of Ag v , introduced in Eq. 

A first way to solve the system (HH1) - ([ST|) has been described in Ref. [l8| and uses evolution equations for h rr 
and //, from which other scalar fields are deduced through the gauge equations (f53"l) - (|55)) as solutions of the angular 
Laplace operator, with radial derivatives as sources. However, this method has the great disadvantage of requiring the 
computation of two radial derivatives to get h lJ , when the source S 13 already contains second-order radial derivatives 
of h %3 . This fourth-order derivation introduces a great amount of numerical noise, which has been observed to rapidly 
spoil the numerical integration. An alternative way is to evolve two other scalar fields and then to integrate (or solve 
PDEs coming from) the Dirac gauge condition to obtain the others. Unfortunately, this is not possible using only the 
six scalar fields (1521) . but one can devise the following procedure in a similar spirit. 

Any rank 2 symmetric tensor T^can be split into two pieces: 



T iJ = [LVj + T 13 = V l V 3 + V 3 V 1 + T 13 , (56) 

with T>jT 13 — 0. For a given T 13 the divergence of Eq. ([56]) allows for the determination of the vector V % through the 
elliptic PDE 

V k V k V l + V l VjV 3 = VjT %3 \ (57) 

where V 1 is fixed up to isometries of fy, which are set by the choice of boundary conditions. If we now return to 
the case T %3 — h %3 and consider only asymptotically flat spatial metric defined on M 3 — no holes — the Dirac gauge 
condition (|50[) is equivalent to having V 1 — 0, since there are no Euclidean symmetries vanishing at infinity. If one 
similarly seeks three scalar fields (^4, B, C) such that: 

A = B = C = f i3 ' = 0, (58) 
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one can check that a solution is: 



dX 








dr 


r 






dW 


A 6v W 


_ I? _ 




dr 


2r 


r 


17' 


dh rr 


2,h rr 
H h 


-2A ev 


fdW 


dr 


r 


\ dr 



B = ^-^^-■±-'—, (60) 
0=^ + ^ + 2^^ + ^). (61) 

In the present case where the trace (or the determinant) is given, B and C are actually coupled and it is sufficient to 
consider: 

B = Y,B em Y lm , with 

to recover B and C using the trace. A nice property of A and B is that, when expressed in terms of these potentials 
related to h lJ , the tensor Poisson equation, with F tJ being a symmetric-tensor representing a source: 

Ah lJ = F lJ (63) 

has a rather simple form. Namely, if we define F A and F B as the scalar potentials similar to A and B, but deduced 
from F 13 , a consequence of Eq. (fBU)) is: 

AA = P 4 , 

AS = F S , (64) 



cr 2 9 1 - 

A := — + -— + — A 6ip and A eip Y em := -£(£ - l)Y tm . (65) 
or^ r dr r z 



with 



Obviously, a very similar property holds for the wave equation (|49| . Therefore, a way of solving numerically the 
constrained system of Eqs. (|49|) -(|51 | . by making use of the potentials A and B, is the following. With the source S 1 ^ 
and /i y known at the initial hypersurface, it is possible to deduce the potentials S A and S B of the source and thus 
to advance the potentials A and B of h lJ to next time-step through the evolution equations 

O 7 \ ^ „ B 



^2- A J B = 5 - ( 66 ) 

Then the six scalar fields (|5"2"|) can be computed by solving the PDE system formed by the following five elliptic 
equations: the definitions of A and B, i.e. Eqs. (|5T))) and (|6"2"|) . together with the Dirac gauge conditions (l5"3"l) - (fB"5"| plus 
the trace- free condition (|5ip — used to get To. All the components of can be finally recovered by taking angular 
derivatives of the scalar fields defined in Eqs. (|B"2"|) . With this algorithm, only two scalar potentials, A and B, are 
evolved in time. The whole tensor is deduced from these potentials and the gauge and trace conditions. Note that, 
when decomposing all the scalar fields onto a spherical harmonics function basis, the elliptic system of five PDEs 
described above reduces to a system of coupled ordinary differential equations in the radial coordinate r. 

With either of these approaches (the one described here or that presented in Ref. [IS]) it is possible to evolve two 
scalar potentials using hyperbolic wave-like operators and recover the symmetric tensor h lJ through an elliptic system 
of PDEs obtained from the gauge conditions. A numerical implementation of these techniques being beyond the scope 
of the present article, we have here only exhibited both algorithms in order to show that it is, in principle, possible to 
build-up the whole conformal metric from the gauge conditions, while being consistent with the evolution equations. 
This might inversely be linked toward the property of the Dirac gauge system being preserved by the 3+1 evolution 
system. Future numerical developments in these directions shall certainly bring better insight into the problem. 
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VI. DISCUSSION. 

All evolution formalisms for the resolution of Einstein equations as an initial value boundary problem exploit the 
intrinsic hyperbolicity of Eqs. |T]), although the associated evolution systems are not necessarily hyperbolic from the 
PDE theory point of view [l| . In the present case of the FCF formalism [l8[ , Einstein equations result in a coupled 
elliptic-hyperbolic PDE system. The hyperbolic part PDE evolution system consists of the reduced system, governing 
the evolution of the gravitational degrees of freedom, whereas the elliptic part is formed by the constrained system 
and part of the gauge system (maximal slicing equation). In fact, in the context of the algorithms presented in section 
IVl the elliptic Dirac system, Eqs. ([5*5)) - ([55|) . can be actually seen as a part of the PDE evolution system. In summary, 
the evolution PDE system is formed by the reduced, constraint, and gauge systems, whereas the the fulfillment of 
the subsidiary system, represented by Eq. (91) in Ref. [l8| for /3 l , can be used as a control test of the scheme along 
the evolution. We have carried out a first analysis of the mathematical structure of the PDE evolution system paying 
particular attention to the equations (JTDJ) governing the evolution for the deviation h lJ of the conformal metric from 
the flat fiducial one /y, i.e. — 7 y — / u . Dirac's gauge plays an important role in getting a well defined hyperbolic 
structure. This elliptic gauge is close in spirit and properties to other gauges employed in the literature, like the 
spatial harmonic gauge in [171 ]. the minimal distortion introduced by York & Smarr, the new minimal distortion 
gauge introduced by Jantzen & York, or the numerically motivated pseudo-minimal distortion gau ge b y Nakamura, 
approximate minimal distortion by Shibata or the Gamma freezing (cf. Sees. 9.3. and 9.4 in Ref. [2 a] for a review 
of them). In particular, all of them can be written as elliptic equations on the shift vector /3\ The Dirac gauge fixes 
spatial coordinates in the evolution (including on the initial data, as the spatial harmonic gauge does) up to boundary 
conditions. For boundary conditions (enforced when solving the elliptic PDE for f] 1 ) such that the evolution vector is 
timclike, the Dirac gauge provides a sufficient condition for the strong hyperbolicity of Eq. (TTU]). Moreover, using this 
gauge it is possible to derive a flux vector in terms of which the first-order system of equations, equivalent to (TIT)]) , has 
the structure of a hyperbolic system of conservation laws (with sources). Likewise, the analysis of the characteristics 
sheds light on the prescription of inner boundary conditions on a spacelike inner cylinder, when employing an excision 
approach to black hole evolutions. More generally, maximal and Dirac gauges can be relaxed to admit more general 
gauges, while preserving the hyperbolic properties of the system but possibly complicating the structure of the sources. 

Having said this, it is clear that further analysis is necessary. First, particular attention should be payed to the source 
terms in equation (TIB")) . They can introduce, in the so-called stiff case, new characteristic time scales (relaxation times 
in the language of fluid dynamics) which may be much smaller than the CFL (Courant-Fricdrichs-Lewy) numerical 
time step (see, e.g., [43l. El 14a) ) . In particular, authors in reference [43| have studied general hyperbolic systems with 
supercharacteristic relaxations, and they shown in which conditions a source term can be damping or, on the contrary, 
enforces growth of instabilities. Looking, in our case, at the quantity Rl J (Eq. I|17p). one can notice the presence of 
quadratic terms in the w% ; it suggests that huge spatial gradients of hp can introduce some degree of stiffness in 
the source terms. Second, nothing has been said about the possible outer boundary conditions to be prescribed when 
studying the initial boundary value problem with an outer timelike cylinder. Certainly, in this case the well-posedness 
analysis is more complicate. However, thanks the enforcement of the constraint along the evolution, there is no need 
of devising specific constraint preserving boundary conditions, and Sommerfeld-like conditions as in (4ll . |42| can be 
straightforwardly employed. Third, nothing has been said about the elliptic part and its coupling with the hyperbolic 
subsystem. On the one hand, this coupling is crucial in the overall well-posedness of the problem, as clearly illustrated 
in the inner boundary conditions issue, where inner boundary conditions on the elliptic part determine the ingoing or 
outgoing nature of the characteristics in the hyperbolic part. On the other hand, the analysis of the elliptic system 
by itself represents an outstanding challenge. This is illustrated by the XCTS elliptic system [26|, (2tJ referred to in 
Section [TBI very closely related to the FCF elliptic subsystem. We note that, in this case, no results on existence are 
available and very little is known on uniqueness, where recent numerical [28;!, 29] and analytical works [13, [HJ works 
point toward the essential non-uniqueness of the system (related to a wrong sign in the differential operator of the 
maximal slicing equation). Fourth, nothing has been said about consequences on well-posedness of coupling matter 
equations to the gravitational degrees of freedom. 

Although our analysis is far from being exhaustive, it has the advantage of giving some clues about which numerical 
strategies are the most convenient in order to solve Einstein equations in the FCF formalism. In this sense, we have 
attempted to obtain some limited but concrete results, rather than remained frozen by the "non-attainability" of 
complete and fully rigorous results. 
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